DP还能干这个？深度剖析DeePMD在溶液反应中的应用（三）
第三回 DPRc渡准劫，gwTP定乾坤
上回书说到，彼时的DFTB2 QM/MM+DPRc[1-2]（下称DPRc），自由能曲线（FES）与PBE0/6-31G* QM/MM（下称PBE0）相比，虽已非常接近，但仍有些许差距。罗格斯大学的Darrin M. York课题组召集众研究者，要还DPRc一个更加精准的FES。传说，研究者讨论时，曾言到，可以把同步学习（concurrent learning）流程中生成的多个DPRc模型与Pengfei Li等人提出的加权热力学扰动（wTP）方法[3]结合起来，获得更加准确的自由能。
据此，研究者基于wTP方法[3]，提出了广义加权热力学扰动（generalized weighted thermodynamic perturbation, gwTP）方法，可以根据多个势能的轨迹，推测出另一个势能的自由能曲线，相关论文以Multireference Generalization of the Weighted Thermodynamic Perturbation Method为题在J. Phys. Chem. A发表[4]。
wTP vs gwTP
研究者首先选取了一个最简单的体系来比较wTP方法和gwTP方法，利用DP-GEN[5]生成了四个DPRc模型，每个模型分别训练了一条FES上的不同区域。图1左侧红线即为DPRc模型的FES，而右侧红线则为wTP对该模型校正后的结果，下方则显示了重加权熵（reweighting entropy, RE），可以理解为权重的“平坦度”。
图1 对某磷酸酯反应使用DPRc产生的自由能曲线。黑线为MABR从PBE0轨迹生成的参考值，左侧红线是MBAR从DPRc轨迹产生的自由能曲线，右侧红线是用wTP从DPRc轨迹估计PBE0产生的自由能曲线，每张图的下方为重加权熵（RE）
可以发现，图1右侧，wTP（红线）和参考值（黑线）在训练区域接近，但在未训练区域相差较远。当红线和黑线十分接近时，对应的重加权熵应当高于0.6。而如果重加权熵降低到0.3左右，那么红线和黑线往往差距甚远。这意味着重加权熵可以验证wTP生成的FES的可靠性，wTP的可靠性依赖于模型在某一构象的准确性。
gwTP把4个模型产生的轨迹结合起来，即“聚合采样”，可以得到图2的曲线，gwTP（红线）和PBE0（黑线）的曲线十分接近，重加权熵也很高，表明gwTP无需每个模型在每一构象都是准确的，某处只要有任一模型准确，足以产生准确的FES。
图2 gwTP产生的自由能曲线。黑线为PBE0轨迹的MBAR参考值，左侧红线仅用DPRc轨迹估计PBE0，右侧同时用DPRc轨迹和PBE0轨迹估计PBE0，下方为重加权熵
gwTP的这一特性，恰恰与DP-GEN中利用同步学习工作流产生的4个模型十分契合。DP-GEN产生的4个模型未必在某一构型完全准确，可能具有些许差别，可以通过模型偏差辅助判断相似度。而gwTP可以同时利用4个不同的模型，聚合采样，产生更加准确的FES。
非酶磷酸酯交换反应
利用gwTP方法，研究者接续上一回的研究，利用上一回使用DP-GEN获得的4个模型[2]，继续探索水溶液中RNA的非酶磷酸酯交换反应。图3左侧显示了100 ps/window的PBE0的FES（黑）和100 ps（红）、1.2 ns（绿）的DPRc FES，其中黑线和红线是上一回的数据，可以看出DPRc相比PBE0仍有差距，并且100 ps/window的DPRc模拟并未达到平衡状态。
图3右侧则是使用gwTP的结果，红线利用100 ps/window的DPRc轨迹预测PBE0的FES，绿线利用1.2 ns/window的DPRc轨迹预测PBE0的FES，蓝线则把PBE0的轨迹和1.2 ns/window一起用来预测FES。可以看出，在gwTP的帮助下，红线较好地复现了黑线。
图3 左侧为MBAR生成的自由能曲线，右侧为gwTP生成的自由能曲线，下方为重加权熵。黑线、红线、绿线分别为100 ps/window PBE0、100 ps/window DPRc和1.2 ns/window DPRc
然而，右侧的红线与绿线也有1 kcal/mol的差距，而它们的重加权熵（图1下方）是相似的，保持在高位，表明它们都是可靠的。因此，可以推断出，100 ps/window的PBE0同样没有达到平衡状态，即黑线的参考值是不可靠的，导致100ps PBE0 FES和DPRc FES产生较大差距。
那么，平衡状态下收敛的FES应当为何？如图4所示，进一步分析发现，DPRc模拟的前600 ps没有达到平衡状态，而600 ps后达到了平衡状态，因此对后600 ps的轨迹使用gwTP方法（红线），取得了比100 ps PBE0（黑线）更加可靠的FES。100 ps/window PBE0得到的能垒和反应能分别是19.6 kcal/mol和1.5 kcal/mol，而把1.2 ns中后600 ns/window的DPRc+gwTP作为最终参考的FES，能垒和反应能分别是21.1 kcal/mol和3.6 kcal/mol。1.2 ns/window的PBE0的成本高昂，无法接受，而相比之下，使用gwTP和DPRc方法，可以以1/86的时间成本取得平衡状态下的FES。
图4
gwTP方法已集成至ndfes软件包中，可以从https://gitlab.com/RutgersLBSR/fe-toolkit/ 下载。
欲知后事如何，且听下回分解。
参考文献
[1] Zeng, J; Giese, T.; Solen, E.; York, D. J. Chem. Theory Comput., 2021, 17, 6993–7009.
[2] Giese, T.; Zeng, J; Solen, E.; York, D. J. Chem. Theory Comput., 2022, 18, 4304–4317.
[3] Li, P.; Jia, X.; Pan, X.; Shao, Y.; Mei, Y., J. Chem. Theory Comput. 2018, 14, 5583–5596.
[4] Giese, T.; Zeng, J; York, D. J. Phys. Chem. A 2022.
[5] Zhang, Y.; Wang H.; Chen W.; Zeng J.; Zhang L.; Wang H.; E W.; Comput. Phys. Commun. 2020, 253, 107206.
推荐阅读